########################################################################################
########################################################################################
####Appendix: Non-state Atrocities in Capital Cities - ZINB Global Robustness Models####
########################################################################################
########################################################################################
library(MASS) 
library(pscl)
library(foreign)
library(stargazer)
library(mvtnorm)
library(plyr)
library(glmmADMB)
library(MCMCglmm)
library(survival)

# Set working library
setwd("~/Data/Global Analysis/")
# Read in main data
main.data <- read.dta("full.grid.upd.dta")

###################################################################
###Model AI-I: Bayesian random effects at country and gid levels###
###################################################################
###Subset the data without NAs
main.data.ai.i<-na.omit(main.data[,c("incidentnonstatefull", "ccode", "lag_Capital", "loglagttime", "lagcivconflagtemp","loglagross_oil_prod", "loglagttime","loglagppp", "lagincidentsumfull", "lagurban", "loglagbdist1", "gid", "loglagcellarea","lagp_polity2", "loglagpop", "year95", "year96","year97","year98","year99","year00","year01","year02","year03","year04","year05","year06","year07","year08","year09", "year")])
####Create the prior
a.V <- matrix(c(100,0,0,100), ncol=2, byrow = TRUE)
zi.prior <-  list(R = list(V = diag(2), n = 1.002, fix = 2),
                  G = list(G1 = list(V = diag(2), n = 1.002,
                                     alpha.mu = c(0,0),
                                     alpha.V = a.V), 
                           G2 = list(V = diag(2), n = 1.002,
                                     alpha.mu = c(0,0),
                                     alpha.V = a.V)))
###The model
at.zinb.ai.i <- MCMCglmm(incidentnonstatefull ~ trait-1 + at.level(trait, 1):lag_Capital + 
                           at.level(trait, 1):year97+at.level(trait, 1):year98+at.level(trait, 1):year99+at.level(trait, 1):year00+
                           at.level(trait, 1):year01+at.level(trait, 1):year02+at.level(trait, 1):year03+at.level(trait, 1):year04+
                           at.level(trait, 1):year05+at.level(trait, 1):year06+at.level(trait, 1):year07+at.level(trait, 1):year08+
                           at.level(trait, 1):year09 +    
                           trait:lagcivconflagtemp + trait:lagincidentsumfull + trait:lagurban + trait:loglagross_oil_prod +
                           trait:loglagppp + trait:loglagbdist1+
                           trait:lagp_polity2 + trait:loglagpop + trait:loglagttime + trait:loglagcellarea, 
                         data = main.data.ai.i, family = "zipoisson",
                         random = ~us(at.level(trait,1) +
                                        at.level(trait,1):ccode):gid + 
                           us(at.level(trait,2) +
                                at.level(trait,2):ccode):gid,
                         rcov = ~ idh(trait):units,
                         prior = zi.prior,
                         nitt = 100000, burnin = 90000, thin = 50, ###Run 100K iterations, keep last 10K with a thinning of 50
                         verbose = TRUE, pr = TRUE, pl = TRUE)
summary(at.zinb.ai.i)

############################################
###Model AI-II: with military expenditure###
############################################
at.zinb.ai.ii <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + loglagmilex + factor(year)
                          |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + loglagmilex, 
                          dist = "negbin", data = main.data)
summary(at.zinb.ai.ii)
AIC(at.zinb.ai.ii)

############################################################
###Model AI-III: with mountains (data only for 2001-2009)###
############################################################
at.zinb.ai.iii <-zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + lagmnt + factor(year)
                          |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + lagmnt, 
                          dist = "negbin", data = main.data)
summary(at.zinb.ai.iii)
AIC(at.zinb.ai.iii)

#########################################
###Model AI-IV: Cluster SEs by country###
#########################################
at.zinb.ai.iv <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                            year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08 + cluster(ccode)
                          |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                          dist = "negbin", data = main.data)
summary(at.zinb.ai.iv)
AIC(at.zinb.ai.iv)

##########################################
###Model AI-V: Cluster SEs by grid cell###
##########################################
at.zinb.ai.v <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                           year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08 + cluster(gid)
                         |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                         dist = "negbin", data = main.data)
summary(at.zinb.ai.v)
AIC(at.zinb.ai.v)

########################################
###Model AI-VI: Country fixed effects###
########################################
at.zinb.ai.vi <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                            year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08 + factor(ccode)
                          |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                          dist = "negbin", data = main.data)
summary(at.zinb.ai.vi)
AIC(at.zinb.ai.vi)

########################################
###Model AI-VII: Remove US grid cells###
########################################
###Subset the data without US
main.data.ns<-subset(main.data,(main.data$ccode!=2))
###The model
at.zinb.ai.vii <-zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                            year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                          |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                          dist = "negbin", data = main.data.ns)
summary(at.zinb.ai.vii)
AIC(at.zinb.ai.vii)

#############################################
###Model AI-VIII: With distance to capital###
#############################################
at.zinb.ai.viii <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcapdist + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                            + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                              year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                            |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                            + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                            dist = "negbin", data = main.data)
summary(at.zinb.ai.viii)
AIC(at.zinb.ai.viii)

###############################################
###Model AI-IX: District year level analysis###
###############################################
#Read in data by district
dis.dat <- read.dta("full.prov.dta")
#Run the model
at.zinb.ai.ix <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                            year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                          |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                          + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                          dist = "negbin", data = dis.dat)
summary(at.zinb.ai.ix)
AIC(at.zinb.ai.ix)

######################################################
###Model AI-X: Use UCDP GED one-sided violence data###
######################################################
at.zinb.ai.x <- zeroinfl(ucdp_inc ~ lag_Capital + lagcivconflagtemp + lagucdp_inc + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                           year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                         |lagcivconflagtemp + lagucdp_inc + lagurban + loglagross_oil_prod
                         + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                         dist = "negbin", data = main.data)
summary(at.zinb.ai.x)
AIC(at.zinb.ai.x)

############################################################
###Model AI-XI: Freedom of press indicator in both phases###
############################################################
at.zinb.ai.xi <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + 
                        lagfred_press + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                      |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + lagfred_press + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data)
summary(at.zinb.ai.xi)
AIC(at.zinb.ai.xi)

#############################################
###Model AI-XII: Outlier countries removed###
#############################################
###Remove outliers
#Summary atrocities by country
p <- summaryBy(incidentnonstatefull ~ ccode, data=main.data, keep.names = T, FUN=c("sum"))
#See which countries are in the 90th percentile
p1 <- subset(p, incidentnonstatefull>=quantile(p$incidentnonstatefull, 0.9))
p1
#Subset only countries below the 90th percentile (this removes the 18 most atrocity prone countries)
main.data.r <- subset(main.data, (ccode!=100 & ccode!=365 & ccode!=475 & ccode!=490 & ccode!=500 & ccode!=516 & ccode!=517 & ccode!=520 & ccode!=615 & ccode!=625 & ccode!=645 & ccode!=666 & ccode!=700 & ccode!=750 & ccode!=770 & ccode!=771 & ccode!=780 & ccode!=850))
###Run the model
at.zinb.ai.xii <- zeroinfl(incidentnonstatefull ~ lag_Capital + lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea + year96 + year97+
                        year98 + year99 + year00 + year01 + year02 + year03 + year04 + year05 + year06 + year07 + year08
                      |lagcivconflagtemp + lagincidentsumfull + lagurban + loglagross_oil_prod
                      + loglagppp + loglagbdist1 + lagp_polity2 + loglagpop + loglagttime + loglagcellarea, 
                      dist = "negbin", data = main.data.r)
summary(at.zinb.ai.xii)
AIC(at.zinb.ai.xii)

###Export to LaTex
#First table
stargazer(at.zinb.ai.ii, at.zinb.ai.iii, at.zinb.ai.iv, at.zinb.ai.v, at.zinb.ai.vi)
stargazer(at.zinb.ai.ii, at.zinb.ai.iii, at.zinb.ai.iv, at.zinb.ai.v, at.zinb.ai.vi, zero.component=TRUE)
AIC(at.zinb.ai.ii, at.zinb.ai.iii, at.zinb.ai.iv, at.zinb.ai.v, at.zinb.ai.vi)
#Second table
stargazer(at.zinb.ai.vii, at.zinb.ai.viii, at.zinb.ai.ix, at.zinb.ai.x, at.zinb.ai.xi, at.zinb.ai.xii)
stargazer(at.zinb.ai.vii, at.zinb.ai.viii, at.zinb.ai.ix, at.zinb.ai.x, at.zinb.ai.xi, at.zinb.ai.xii, zero.component=TRUE)
AIC(at.zinb.ai.vii, at.zinb.ai.viii, at.zinb.ai.ix, at.zinb.ai.x, at.zinb.ai.xi, at.zinb.ai.xii)


